# -------------------------------------------------------------------
# Purpose: Creates Table B43
# Author:  Max Posch, 25/07/2025
# Usage:   Source this script to generate the table.
# -------------------------------------------------------------------
# Check that required paths exist
stopifnot(dir.exists(pdataanalysis))
stopifnot(dir.exists(poutputappendix))


# Load data
load(file.path(pdataanalysis, "countyLevel19001940.RData"))
d <- copy(countyLevel19001940)[year==1940]


## Winsorize and normalize
vars <- str_subset(names(d), "^sum[:graph:]+pc_f$")
d[, paste0(vars, "_w") := lapply(.SD, function(x) winsorize(x, probs = c(.01, .99))), .SDcols = vars]

vars <- str_subset(names(d), "entropy|sum_n_namelast_mp_adjp|mean_yschool_adjp")
d[, paste0(vars, "_w") := lapply(.SD, function(x) winsorize(x, probs = c(.01, .99))), .SDcols = vars]

vars <- str_subset(names(d), "entropy|sum_n_namelast_mp_adjp|mean_yschool_adjp")
vars <- str_subset(vars, "_w$")
d[, paste0(vars, "s") := lapply(.SD, scale), .SDcols = vars]


# Regressions
o <- list()
o <- append(o, list(feols(sum_patents_pc_f_w ~ mean_yschool_adjp_ws + sum_n_namelast_mp_adjp_ws | statefip, d)))
o <- append(o, list(feols(sum_patents_pc_f_w ~ entropy_namelast_mp_adjp_ws + mean_yschool_adjp_ws + sum_n_namelast_mp_adjp_ws | statefip, d)))
o <- append(o, list(feols(sum_patents_pc_f_w ~ (entropy_namelast_mp_adjp_ws + mean_yschool_adjp_ws)^2 + sum_n_namelast_mp_adjp_ws | statefip, d)))
o <- append(o, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ mean_yschool_adjp_ws + sum_n_namelast_mp_adjp_ws | statefip, d)))
o <- append(o, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ entropy_namelast_mp_adjp_ws + mean_yschool_adjp_ws + sum_n_namelast_mp_adjp_ws | statefip, d)))
o <- append(o, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ (entropy_namelast_mp_adjp_ws + mean_yschool_adjp_ws)^2 + sum_n_namelast_mp_adjp_ws | statefip, d)))

r <- list()
r <- append(r, list(feols(sum_patents_pc_f_w ~ mean_yschool_adjp_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | statefip, d)))
r <- append(r, list(feols(sum_patents_pc_f_w ~ iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + mean_yschool_adjp_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | statefip, d)))
r <- append(r, list(feols(sum_patents_pc_f_w ~ (iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + mean_yschool_adjp_ws)^2 + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | statefip, d)))
r <- append(r, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ mean_yschool_adjp_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | statefip, d)))
r <- append(r, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + mean_yschool_adjp_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | statefip, d)))
r <- append(r, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ (iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + mean_yschool_adjp_ws)^2 + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | statefip, d)))

i <- list()
i <- append(i, list(feols(sum_patents_pc_f_w ~ mean_yschool_adjp_ws | statefip | sum_n_namelast_mp_adjp_ws ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws, d)))
i <- append(i, list(feols(sum_patents_pc_f_w ~ mean_yschool_adjp_ws | statefip | entropy_namelast_mp_adjp_ws + sum_n_namelast_mp_adjp_ws ~ iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws, d)))
i <- append(i, list(feols(sum_patents_pc_f_w ~ mean_yschool_adjp_ws | statefip | entropy_namelast_mp_adjp_ws + sum_n_namelast_mp_adjp_ws + entropy_namelast_mp_adjp_ws:mean_yschool_adjp_ws ~ iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws + iv_lo_entropy_namelast_mp_adjp_fe_immig_ws:mean_yschool_adjp_ws, d)))
i <- append(i, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ mean_yschool_adjp_ws | statefip | sum_n_namelast_mp_adjp_ws ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws, d)))
i <- append(i, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ mean_yschool_adjp_ws | statefip | entropy_namelast_mp_adjp_ws + sum_n_namelast_mp_adjp_ws ~ iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws, d)))
i <- append(i, list(feols(sum_break_p80_rrfsim05_pc_f_w ~ mean_yschool_adjp_ws | statefip | entropy_namelast_mp_adjp_ws + sum_n_namelast_mp_adjp_ws + entropy_namelast_mp_adjp_ws:mean_yschool_adjp_ws ~ iv_lo_entropy_namelast_mp_adjp_fe_immig_ws + iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws + iv_lo_entropy_namelast_mp_adjp_fe_immig_ws:mean_yschool_adjp_ws, d)))


# F-statistics
dstata <- d[, .(
  x1 = entropy_namelast_mp_adjp_ws, z1 = iv_lo_entropy_namelast_mp_adjp_fe_immig_ws,
  x2 = sum_n_namelast_mp_adjp_ws, z2 = iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws,
  y1 = sum_patents_pc_f_w, y2 = sum_break_p80_rrfsim05_pc_f_w,
  w = mean_yschool_adjp_ws,
  statefip_f, year_f
)]
dstata[, x1_ws := x1 * w]
dstata[, z1_ws := z1 * w]
commands <- list(
  'ivreghdfe y1 w (x2 = z2), absorb(statefip) cluster(statefip_f) first ffirst savefirst
   gen swf1 = .
   replace swf1 = round(e(first)["SWF",1])
   keep swf*
   keep if _n == 1',
    'ivreghdfe y1 w (x1 x2 = z1 z2), absorb(statefip) cluster(statefip_f) first ffirst savefirst
   gen swf1 = .
   gen swf2 = .
   replace swf1 = round(e(first)["SWF",1])
   replace swf2 = round(e(first)["SWF",2])
   keep swf*
   keep if _n == 1',
    'ivreghdfe y1 w (x1 x2 x1_ws = z1 z2 z1_ws), absorb(statefip) cluster(statefip_f) first ffirst savefirst
   gen swf1 = .
   gen swf2 = .
   gen swf3 = .
   replace swf1 = round(e(first)["SWF",1])
   replace swf2 = round(e(first)["SWF",2])
   replace swf3 = round(e(first)["SWF",3])
   keep swf*
   keep if _n == 1'
)
swf_results <- as.character(unlist(lapply(commands, get_fstat_from_stata, data.in = dstata)))
swfstat <- rep(swf_results, times = 2)


# Create table
y1_mean <- round(mean(d$sum_patents_pc_f_w), 2)
y1_sd <- round(sd(d$sum_patents_pc_f_w), 2)
y2_mean <- round(mean(d$sum_break_p80_rrfsim05_pc_f_w), 2)
y2_sd <- round(sd(d$sum_break_p80_rrfsim05_pc_f_w), 2)
y1 <- paste0("\\makecell{Patents \\\\ per 1,000 people \\\\ (mean = ", y1_mean, ", sd = ", y1_sd, ")}")
y2 <- paste0("\\makecell{Breakthrough patents \\\\ per 1,000 people \\\\ (mean = ", y2_mean, ", sd = ", y2_sd, ")}")

setFixest_dict(
  c(
    entropy_namelast_mp_adjp_ws = "Surname diversity",
    iv_lo_entropy_namelast_mp_adjp_fe_immig_ws = "Predicted surname diversity",
    sum_n_namelast_mp_adjp_ws = "Population",
    iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws = "Predicted population",
    mean_yschool_adjp_ws = "Average years of schooling",
    sum_patents_pc_f_w = y1, sum_break_p80_rrfsim05_pc_f_w = y2,
    year = "Period", statefip = "State", gisjoin = "County"
  )
)

tablename <- file.path(poutputappendix, "tableB43.tex")
etable(o,
  cluster = ~statefip,
  fitstat = ~n,
  order = c("Surname diversity$", "^Average years of schooling", "Surname", "Population"),
  digits = "r3", digits.stats = "r3",
  file = tablename, replace = TRUE,
  style.tex = style.tex("aer"), tex = TRUE
)
edit_table_content_fixed(tablename, "Period $\\times $ County", "County-specific linear trends")
add_table_row(tablename, "\\midrule", "\\multicolumn{2}{l}{\\textit{Panel A: Least-squares estimates}} &  \\multicolumn{5}{c}{}\\\\ \\cmidrule(lr){1-7}")
add_table_row(tablename, "mean =", "\\cmidrule(lr){2-4}  \\cmidrule(lr){5-7}")
move_table_row(tablename, "Observations", "bottomrule")
add_table_row(tablename, "    \\\\", c("\\multicolumn{2}{l}{\\textit{Panel B: Reduced-form estimates}} &  \\multicolumn{5}{c}{}\\\\", "\\\\", "\\multicolumn{2}{l}{\\textit{Panel C: Instrumental-variable estimates}} &  \\multicolumn{5}{c}{}\\\\", "\\\\", "\\multicolumn{2}{l}{\\textit{Panel D: First-stage estimates}} &  \\multicolumn{5}{c}{}\\\\"))

temptable <- file.path(poutputappendix, "temp.tex")
etable(r,
  cluster = ~statefip,
  fitstat = ~n,
  order = c("urname diversity$", "^Average years of schooling", "urname", "opulation"),
  digits = "r3", digits.stats = "r3",
  file = temptable, replace = TRUE,
  style.tex = style.tex("aer"), tex = TRUE
)
estimates_rows <- get_estimates_rows(temptable)
add_table_row(tablename, "Panel B", c("\\cmidrule(lr){1-7}", estimates_rows))

temptable <- file.path(poutputappendix, "temp.tex")
etable(i,
  cluster = ~statefip,
  fitstat = ~n,
  order = c("Surname diversity$", "^Average years of schooling", "Surname", "Population"),
  digits = "r3", digits.stats = "r3",
  extralines = list("Sanderson-Windmeijer \\textit{F}-stat" = swfstat),
  file = temptable, replace = TRUE,
  style.tex = style.tex("aer"), tex = TRUE
)
estimates_rows <- get_estimates_rows(temptable)
add_table_row(tablename, "Panel C", c("\\cmidrule(lr){1-7}", estimates_rows))
add_table_row(tablename, "Sanderson", "\\\\", "before")
edit_table_content_fixed(tablename, "entropy\\_namelast\\_mp\\_adjp\\_ws:mean\\_yschool\\_adjp\\_ws", "Surname diversity $\\times$ Average years of schooling")

temptable <- file.path(poutputappendix, "temp.tex")
etable(i,
  stage = 1,
  cluster = ~statefip,
  order = c("urname diversity$", "^Average years of schooling", "urname", "opulation"),
  fitstat = ~n, digits = "r3", digits.stats = "r3",
  file = temptable, replace = TRUE,
  style.tex = style.tex("aer"), tex = TRUE
)
estimates_rows <- get_estimates_rows(temptable)
estimates_rows1 <- collapse_stage1(estimates_rows, c(NA, 3, 5,  2, 4, 6))
add_table_row(tablename, "Panel D", estimates_rows1)
add_table_row(tablename, "Panel D", c("\\cmidrule(lr){1-7}", "& \\multicolumn{3}{c}{Surname diversity} & \\multicolumn{3}{c}{Population}\\\\", "\\cmidrule(lr){2-4}  \\cmidrule(lr){5-7}"))
estimates_rows2 <- collapse_stage1(estimates_rows, c(NA, NA, 7, NA, NA, NA))
add_table_row(tablename, "State fixed", c("& \\multicolumn{3}{c}{\\makecell{Surname diversity $\\times$ Average years of schooling}} & \\multicolumn{3}{c}{}\\\\", "\\cmidrule(lr){2-4}", estimates_rows2, "\\\\"), "before")
file.remove(temptable)

cat("Table B43 saved to:", tablename, "\n")